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Abstract. We present a numerical method and computer 
code to calculate the radiative transfer and excitation 
of molecular lines. Formulating the Monte Carlo method 
from the viewpoint of cells rather than photons allows us 
to separate local and external contributions to the radia- 
tion field. This separation is critical to accurate and fast 
performance at high optical depths (r ^ 100). The ran- 
dom nature of the Monte Carlo method serves to verify the 
independence of the solution to the angular, spatial, and 
frequency sampling of the radiation field. These features 
allow use of our method in a wide variety of astrophysi- 
cal problems without specific adaptations: in any axially 
symmetric source model and for all atoms or molecules for 
which collisional rate coefficients are available. Continuum 
emission and absorption by dust is explicitly taken into 
account but scattering is neglected. We illustrate these 
features in calculations of (i) the HCO+ J=l-0 and 3-2 
emission from a fiattened protostellar envelope with infall 
and rotation, (u) the CO, HCO+, CN and HCN emission 
from a protoplanetary disk and (iii) HCN emission from 
a high-mass young stellar object, where infrared pump- 
ing is important. The program can be used for optical 
depths up to 10'^ — lO**, depending on source model. We 
expect this program to be an important tool in analysing 
data from present and future infrared and (sub) millimetre 
telescopes. 

Key words: Line: formation - Radiative transfer - Meth- 
ods: numerical - Stars: formation - ISM: molecules 



1. Introduction 

The dense and cool material in the interstellar medium 
of galaxies plays an important role in the life cycle of 
stars, from the earliest phases of star formation to the 



shells around evolved stars and the gas and dust tori 
around active galactic nuclei. Line emission from atoms 
and molecules, and continuum emission from dust parti- 
cles, at radio, (sub) millimetre and infrared wavelengths 
are indispensable tools in the study of a wide variety of as- 
trophysical problems. This is illustrated by the large num- 
ber of infrared and submillimetre observatories planned 
for the near future, such as the Smithsonian Millime- 
ter Array (SMA), the Atacama Large Millimeter Array 
(ALMA), the Far- Infrared and Submillimetre Space Tele- 
scope (FIRST) and the Stratospheric Observatory for In- 
frared Astronomy (SOFIA). 

An essential step in the interpretation of the data from 
these instruments is the comparison with predicted emis- 
sion from models. This paper presents a numerical method 
to solve the radiative transfer and molecular excitation in 
spherically symmetric and cylindrically symmetric source 
models. At the comparatively low densities of interstellar 
gas, the excitation of many molecules is out of local ther- 
modynamic equilibrium (LTE), and the transfer of line 
(and continuum) radiation plays a significant role in detcr- 
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mining the molecular excitation (Leung and Liszt, 1976; 
p31ack, 2000| ). Geometry thus becomes an important ele- 
ment, and the high spatial resolution of current and future 
instruments often demands that at least two-dimensional 
(axisymmetric) source structures are considered. In the 
implementation of our method discussed in this paper, we 
have limited the source structure to spherical and cylindri- 
cal symmetries. The large and often multidimensional pa- 
rameter space further requires a fast and reliable method, 
which needs to be easily applicable to many different as- 
trophysical problems. 

This need for reliable and fiexible tools calls for the 
use of Monte Carlo techniques, where the directions of 
integration are chose n ran domly. This approach was first 
explored by Bernes ( 197!:) for non-LTE mole cular exci- 
tation; later, Choi et al. ( 1995 ), Juvela (1997) and Park 
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& Hong ( 199§| ) augmented it and expanded it to multi- 
ple dimensions. However, Monte Carlo methods can be 
quite slow, especially at large optical depths (r > 100), 
which has limited their use so far. We will show that this 
problem can be overcome by using a technique inspired 
on Accelerated Lambda Iteration: the local radiation field 
and excitation are solved self-consistently and separated 



from the overall radiative transfer problem (see §3.4). The 
greatest virtue of our code is its ability to deal with a wide 
variety of source models for many atomic and molecular 
species, with or without a dust continuum. Although for 
any individual problem a somewhat more efficient method 
could be constructed (§ 3^), the Monte Carlo approach 
frees the user from having to fine-tune such a method and 
allows the user to focus on the astrophysics of the problem 
at hand. 

This paper does not discuss the influence of radiative 
transfer on the source structure, through the thermal bal- 
ance, ionization and chemistry ( Takahashi ct al., 1983; 
Ceccarelh ct al., 199^ ; ^oty and Neufeld, 1997| ,^ exam- 
ple). However, our code is suited to become part of an 
iterative scheme to obtain self-consistent solutions for the 
source structure including radiative transfer and molecu- 
lar excitation. 

Throughout this paper, examples from studies of star 
formation will serve to illustrate the various topics - and 
to stress the link with the analysis of observations. Sect. || 
introduces a simple, spherically symmetric model of the 
core of an interstellar cloud, collapsing to form a star. 
Sect. H then discusses the coupled problem of radiative 
transfer and molecular excitation. It introduces the two 
most commonly used solution methods, and shows that 
these are closely related. This opens the possibility of con- 
structing a hybrid method which combines the benefits of 
both; the implementation of this combined approach in 
our code is deferred to the Appendix. The paper continues 
by exploring the capabilities of our code through a number 
of astrophysically relevant examples, based on extensions 
of the simple one-dimensional model of Sect. ||. A brief 
summary concludes the paper in Sect. ^. 

2. An illustrative, one-dimensional model 

The formation of stars occurs in dense condensations 
within interstellar molecular clouds, which collapse under 
the influence of their own gravity. A widely used theoreti- 
cal description of this process, constructed by Shu ( 1977 ), 
starts with the singular isothermal sphere. 



p(r,t = 0) 



27rG 



(1) 



Here, p is the density as function of radius r and time t, a 
is the isothermal sound speed, and G is the gravitational 
constant. 

At i = collapse starts at the center (r = 0). Af- 
ter a time t, all regions r < at are collapsing, with speed 
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Fig. 1. Density (top left) and velocity (bottom left) 
structure of the spherically-symmetric inside-out collapse 
model of Shu (1977) used to illustrate the radiative trans- 
fer and molecular excitation problem (S|2]). The excitation 
of HCO+ (top right; solid line) ranges from LTE in the 
dense, central regions to sub-thermal in the lower density 
outer regions. Compared to the optically thin excitation 
of H^'^CO"'" (top right; dashed line), line trapping signif- 
icantly influences the HCO+ excitation. The distribution 
of the kinetic temperature is shown with the thick line 
for comparison. The lower right panel shows the emergent 
HCO+ and Hi3C0+ J=4-3 line profiles in a 14" beam for 
a source at 140 pc. The asymmetric profile of the optically 
thick HCO+ 4-3 line is characteristic of infall. 
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v(r,t) increasing from at r = at to free-fall, v oz r 
well within this 'collapse expansion wave' (r <C at). Shu 
( 1977D constructed a solution for the density and veloc- 
ity field of the collapsing core which is self-similar in the 
spatial coordinate x = at. The density follows a power-law 
behaviour as function of radius, with p cx r~^'^ for r <C at, 
p cx just inside r — at, and the undisturbed p cx 
outside r — at (Fig. |l]). 

Many authors have tested this model against obser- 
vations of cloud cores and envelopes around young stellar 
objects (YSOs), e.g. Zho u et a l. ( |1993| ), Choi et al. ( [1995D , 
Ward-Thompson et al. (1996) and Hogerheijde & Sandell 
( [2000D. Especially the spectral-line signature of collapse 
(Fig. Illd) has received much attention as a probe of on- 
going collapse, although this signature is shared by all 
collapse models and is not unique to the particular model 
described here. The exact line shape, however, depends 
quantitatively on the adopted model. 
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Table 1. Molecular data used in this paper 



Molecule No. Levels References 



CO" 
CO' 
HCO+ 

cs 

CN 

HCN 

0-H2CO 



6 
26 
21 
12 
15 
36"^ 
20 



Green fc Thad deus ( [1976| ) 
Schinke et al. ( | l985|) 
Monteiro ( |l98g ) 
Green & Ch apma n ( |l97^ ) 
Black et al. ( [l99l| ) 



Green (1994, priv. comm.)*^ 
Green (|l99l|) 



a Calculation presented in A ppen dix | 
b Calculation presented in §§ [l.l| and I 
c Based on results of Green & Chapman ( 1978 ) for CS. 
d Lev els up to J = 10 in both the 1^2 = and V2 = 1 states. 




See http:/ /www. giss.nasa.gov/data/mcrates 



The interpretation of this signature needs non-LTE ra- 
diative transfer. Both colhsional and radiative processes 
can excite molecules, and for each transition a critical 
density can be defined where the two are of equal im- 
portance. At lower densities radiation dominates, while 
at higher densities collisions drive the level populations 
to thermodynamic equilibrium. The large range of den- 
sities of star forming cores ensures that many molecules 
and transitions will go through the entire range of excita- 
tion conditions, while line emission will have a significant 
impact on the excitation at the intensities and opacities 
expected for typical abundances of many species, not only 
locally but throughout the envelope (Fig. ||c). 

In the following we will use this model to illustrate 
our method of solving the coupled problem of radiative 
transfer and excitation. In particular, we will consider 
emission lines of IICO+ and H^^CO^, which are read- 
ily observed and often used as tracers of dense gas. The 
strong J=1^0,3-^2 and 4^3 lines at millime- 
tre wavelengths have critical densities of 2 x 10^, 4 x 10^, 
and 1 X 10^ cm~'^, using the molecular data in Table |l|. 
We assume an abundance of IIC0"'"/Il2 = 5 x 10~^ and 
an isotopic ratio of 1:65 for H"CO+:HCO+. The sound 
speed of the adopted model is a = 0.24 km s^^, its age 
t = 1 X 10^ yr, and its outer radius 8000 AU. The to- 
tal mass of the model is 0.73 M©. The kinetic tempera- 
ture follows Tkin = 30K(r/1000AU)~° '^, appropriate for 
a centrally heated envelope at a luminosity of ^ 2 L© 
(Adam p et al., 19 8? , e.g.). The turbulent line width of 0.2 



3.1. The coupled problem 

The equation of radiative transport reads, in the notation 
of Rybicki & Lightman ( |l979D , 



dl, 
ds 

or, equivalently, 

dl^ 

dTi, 



(2) 



(3) 



Here, is the intensity at frequency v along a particular 
line of sight parameterized by ds, is the absorption 
coefficient in units cm~^, and ji^ the emission coefficient 
with units erg s~^ cm~^ Hz~^ sr^^. The second form of the 
equation is a useful change of variables, with the source 
function Si, = ju/oii, and the optical depth dT^ = a^ds. 
We consider both molecules and dust particles as sources 
of emission and absorption (j^, = j,y(dust) -I- jy(gas); a^, = 
ay(dust) -I- a,y(gas)), but ignore scattering. Although not 
impossible to include in our code, scattering effects are 
usually negligible at wavelengths longer than mid-infrared. 

When a,y and j^, are known at each position in the 
source, the distribution of the emission on the sky simply 
follows from ray tracing. However, in many cases, and 
ji, will depend on the local mean intensity of the radiation 
field 



J. 



1 

47r 



(4) 



Here, Ji, is the average intensity received from all solid 
angles dQ, and is the solution of Eq. (|^) along each 
direction under consideration. The latter integration ex- 
tends formally to infinity, but in practice only to the edge 
of the source with any incident isotropic radiation field 
like the cosmic microwave background (CMB) as bound- 
ary condition. 

For thermal continuum emission from dust, j;y(dust) 
and ay(dust) are simply given by 



j^(dust) = a,,(dust)B,,(Tdust), 



(5) 



where B^, is the Planck function at the dust temperature 
Tdust, and 



a,y(dust) = Ki^pdust, 



(6) 



where k^, is the dust opacity in cm ^ per unit (dust) mass 
and pdust is the mass density of dust. Our code can use any 



km s ^ is amaller than the ayotcmatic volocitieo except in 
the outermost part (Fig. |^b). 



3. Solving radiative transfer and molecular 
excitation 



et al., 1994 



" Ki, (Osscnkopf and Hcnning, 1994; Pollack 


Draine and Lee, 1984 




Mathis et al., 1977, 



e.g.). 

In the case of emission and absorption in a spectral 
line, a"' (gas) and j"'(gas) are determined by absorption 
and emission between radiatively coupled levels u and I 
with populations (in cm^"^) n.^ and n;. The energy differ- 
ence between levels AE = — Ei corresponds to the 
rest frequency of the transition, = AE/h, where h 
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is Planck's constant. The emission and absorption coef- 
ficients between levels u and / are strongly peaked around 
Uq with a frequency dependence described by a line-profile 
function 



J"'(gas) 
a"' (gas) 



47r 
47r 



nuAui(t){iy), 

{niBiu - nuBui)4'(,y). 



(7) 
(8) 



The Einstein A^i , Biu , and Bui coefficients determine the 
transition probabilities for spontaneous emission, absorp- 
tion, and stimulated emission, respectively, and depend 
on molecule. In most interstellar clouds the line profile is 
dominated by Doppler broadening due to the turbulent 
velocity field 



exp 



(9) 



where the turbulence is assumed to be Gaussian with a full 
width h. In the presence of a systematic velocity field, the 
line profile is angle-dependent and the projection of the 
local velocity vector onto the photon propagation direction 
enters (y — vp). 

Together, collisions and radiation determine the level 
populations through the equation of statistical equilib- 



(10) 



Yl,k>i "^kAki + J2k^i nkiBkiJu + Cki). 

The collision rates Cki depend on the density and the 
coUisional rate coefficients of molecular hydrogen and 
other collision partners, and on temperature through the 
detailed balance of the up- and downward coefficients. 
Eq. (^0|) can be easily solved through matrix inversion 
for each position in the source provided the radiation field 

is known. However, contains contributions by the 
CMB, dust and spectral lines, and since the spectral fine 
term depends on the level populations through Eqs. (^, 
(^ and (||), the problem must be solved iteratively. Start- 
ing with an initial guess for the level populations, J^, is 
calculated, statistical equilibrium is solved and new level 
populations are obtained; through the Monte Carlo in- 
tegration, the new populations yield a new value for J^,, 
after which the populations are updated; etc., until the 
radiation field and the populations have converged on a 
consistent solution. 

When the physical conditions do not vary much over 
the model, an approximate value of J^, can be found 
from the local conditions alone. This idea is the ba- 
sis of the Large Velocity Gradient method, the Sobolev 
method, microturbulence, or the escape probability for- 



malism (3obolev, 196C; Goldreich and Kwan, 1974; Leung 



and Liszt, 1976; de Jong et al., 198C, e.g.). Also, in spe- 



cific geometries, the integration over all solid angles and 
along the full length of the line of sight of Eqs. (||) and 



(^ can be greatly reduced, making the problem tractable. 
This sort of technique has most application in stellar and 
planetary atmospheres; the Eddington approximation is 
an example. 

However, in many astrophysical situations including 
the example of § |^, such simplifications cannot be made, 
and Eqs. (||) and (^ need to be fully solved to get Jy. 
Compared to the relative ease with which statistical equi- 
librium can be solved (Eq. |l^) , obtaining becomes the 
central issue. Direct integration of Eqs. (§) and (^ with, 
e.g., Romberg's method, is infeasible for realistic models, 
but based on a finite set of directions a good approxima- 
tion of Jjy can be obtained. The next two sections describe 
two different methods to choose this set and construct 
in this way. 



3.2. Constructing Ji, and the A-operator 

For computational purposes, source models are divided 
into discrete grid cells, each with constant properties (den- 
sity, temperature, molecular abundance, turbulent line 
width, etc.). It is also assumed that the molecular exci- 
tation can be represented by a single value in each cell, 
which requires instantaneous spatial and velocity mixing 
of the gas. Appropriate source models have small enough 
cells that the assumption of constant excitation is valid. 
The systematic velocity field is the only quantity that is a 
vector field rather than a scalar field, and in our code it is 
allowed to vary in a continuous way within each cell. We 
divide the integration along a ray into subunits within a 
cell to track the variation of the velocity projected on the 
ray. 

Such a gridded source model lends itself easily to the 
construction of a finite set of integration paths to build up 
Ji,. The average radiation field can be thought of as the 
sum of the emission received in cell i from each of the other 
cells j after propagation through the intervening cells and 
weighted with the solid angle subtended by each of these 
cells j as seen from cell i. The combination of radiative 
transfer and statistical equilibrium can be written as 



J, = A[SM)]. 



(11) 



This equation states that the radiation field is given by 
an operator A acting on the source function Sui, which 
depends on the level populations and hence Ji, (Eqs. 0, |[ 
p^ ) . Considering the narrow frequency interval around the 
transition u-l, we have replaced S^, by Sui = [j,yo(dust) -I- 
/ j"'(gas)(i:/]/[Q!i/(,(dust)-|- / a"'(gas)(i:^]. This corresponds 
to the assumption of instanteous redistribution of excita- 
tion mentioned above. In a gridded source model, one can 
think of A as a matrix describing how the radiation field in 
cell i depends on the excitation in all other cells. The ele- 
ments in the matrix then represent the radiative coupling 
between cell pairs. 
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a) Monte Carlo method, taking photon point of view. 
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b) Monte Carlo method, taking cell point of view. 

Fig. 2. (a) In the 'traditional' formulation of the Monte 
Carlo method for solving radiative transfer, the radiation 
field is represented by a certain number of photon pack- 
ages. Each of these packages originates in a random posi- 
tion of the cloud, corresponding to spontaneous emission, 
and travels in a random direction through the cloud until 
it either escapes or is absorbed. To include the CMB field, 
a separate set of packages is included, shown as dashed 
arrows, that originate at the cloud's edge. The packages 
traversing a cell during an iteration give in that cell, 
(b) In our implementation, an equivalent estimate of Ji, 
is found by choosing a certain number of rays which en- 
ter the cell from infinity (or the cloud's edge, using the 
CMB field as a boundary condition) from a random di- 
rection and contribute to the radiation field at a random 
point in the cell's volume. As § |3.4| argues, this formula- 
tion allows separation between the incident radiation field 
and the locally produced radiation field, which accelerates 
convergence in the presence of significant optical depth. 



Eq. ( [Ill ) can be solved iteratively, where an updated 
value of Ji, is obtained by having A operate on the previous 



populations, S 



t 

ul ' 



Since is already known, this only involves matrix mul- 
tiplication, compared to the much more expensive matrix 
inversion required to solve Eq. ( pT] ) . Because of this elegant 
notation, iterative schemes for non-LTE excitation and ra- 
diative transfer are commonly referred to as A-iteration, 
even if no A-operator is ever actually constructed. These 
methods share the use of the same set of rays through- 
out the calculation, in contrast to Monte Carlo methods, 
which use random rays as discussed in § ^3.3| and 3.5, 

Constructing the A-operator in multidimensional 
source models is taxing on computer memory because of 
all the possible connections to keep track of simultane- 
ously. Techniques exist to reduce the number of elements 



(DuUemond and TuroUa, 2000), but these are complex and 
may require some fine-tuning for individual source geome- 
tries. Alternatively, computer memory can be exchanged 
for computing time by solving the problem serially, cal- 
culating the radiation field in each of the cells due to the 
other cells one at a time. 



3.3. The Monte Carlo method 



One way of solving Eq. (y_2|) is to directly sum the contri- 
bution from all other cells to the radiation field in each 
of the individual cells. This corresponds to replacing the 
integral in Eq. by a summation. With a judicially cho- 
sen fixed set of directions or rays, as most A-iteration 
codes do, a good approximation of Ji, can be found in 



J, = K[Sl{J,)]. 



(12) 



this way (Phillips, 1999, e.g.). However, this procedure re- 
quires care, since the necessary angular sampling depends, 
in principle, on characteristics of the excitation solution of 
the problem at hand. 

Since our aim is to construct a method that can be ap- 
plied to many different source models without too much 
fine-tuning, we adopt the Monte Carlo approach to ob- 
tain J^. Analogous to the Monte Carlo method to solve 
the definite integral of a function [see chapter 7 of Press 
et al. (1992) for a discussion of Monte Carlo integration, 
and further references], Eq. (^ can be approximated by 
the summation over a random set of directions. This has 
the advantage that all directions are sampled to sufficient 
detail: if too few directions are included, subsequent re- 
alizations will give different estimates of (see § |3.5| for 
further discussion of this issue). 

Originally ( Bcrncs, 1979| ), the Monte Carlo approach 
was phrased in terms of randomly generated 'photon pack- 
ages', which are followed as they travel through the source 
and which together approximate the radiation field. Fig. |^ 
illustrates that a formulation in terms of randomly chosen 
directions from each cell yields an equivalent estimate of 
Jjy. The only difference is the direction of integration in 
Eq. (|^). Where the former approach follows the photons 
as they propagate through the cells, the latter backtracks 
the paths of the photons incident on each cell. As the 
next section will discuss, this latter approach lends itself 
better to improvements in its convergence characteristics. 
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iteration are large. Large optical depth can be a major 
obstacle to this behaviour: emission passing through an 
opaque cell will rapidly lose all memory of its incident in- 
tensity and quickly tend toward the local source function. 
The distance over which the information about changes 
in excitation can travel is one mean free path per itera- 
tion, so that the required number of iterations grows oc 
characteristic of random walk. This effect makes it hard 
to determine if the process has converged. 



Accelerated or Approximated Lambda Iteration (Ry 



10 

iteration 
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Fig. 3. Evolution of the fractional error in the level popu- 
lations of HCO+ as function of iteration step with ('accel- 
erated'; solid line) and without ('not accelerated'; dashed 
line) separation of local and incident radiation field. For 
the optically thin H^'^CO+ molecule, both methods con- 
verge equally fast. The solid symbols indicate the iteration 
where the first stage of the calculation has converged (see 
§ |3.5| ); after that random noise starts to dominate the frac- 
tional error, which is controlled by the increase in rays per 
cell. The source model is that described in §H. The 'not 
accelerated' HCO"*" formally converged at iteration 18, be- 
cause the difference with iteration 17 became smaller than 
1/30, even though the difference from the real solution ex- 
ceeds that value. This illustrates that acceleration is not 
only computationally convenient, but may also essential 
for a correct solution. 



Treatment of non-isotropic scattering is more complicated 
in this approach, and since scattering is not important at 
the wavelengths of interest here, <; lO/xm, scattering is not 
included in the code. Implementations of the Monte Carlo 
method more appropriate for scattering are available in 
the literature (jWood ct al., 1996^ ; [Wood et al., 1996b| ; 



IWolf ct al. , 19991) 



3.4- Convergence and acceleration 

Besides estimating J^, an important aspect of non-LTE 
radiative transfer is convergence towards the correct solu- 
tion in a reasonable amount of time. Since the solution is 
not a priori known, convergence is often gauged from the 
difference between subsequent iterative solutions. This re- 
lies on the assumption that when and the populations 
are far from their equilibrium solution, corrections in each 



bicki and Hummer, 1991, ALI), circumvents this problem 
by defining an approximate operator A* such that 



J. = iA-A*)[sUj.)] + A*[SuiiJ.)]. 



(13) 



An appropriate choice for A* is one which is easily invert- 
ible and which steers rapidly toward convergence. This 
occurs if Jj/ is dominated by the second term on the right 
hand side of the equation, where A* works on the cur- 
rent source function as opposed to the solution from the 
previous iteration. 



After several attempts (Bcharmer, 1981), Olson, Auer, 
& Buchler (1986) found that a good choice for A* is the 



diagonal, or sometimes tri-diagonal, part of the full oper- 
ator A. This choice for A* describes the radiation field 
generated locally by the material in each cell, and its 
direct neighbours in the case of the tri-diagonal matrix. 
Eq. ( |l3|) then gives as the sum of the field incident on 
each cell due to the previous solution for the excitation 
{(A — A*)[S'^J}, and a self- consistent solution of the local 
excitation and radiation field {A*[S'„;]}. In opaque cells, 
the radiation field is close to the local source function, and 
Eq. ( p^ converges significantly faster than Eq. (^2|); for 
optically thin cells, both formalisms converge equally fast. 

Formulating the Monte Carlo method in terms of ran- 
domly generated photon packages traveling through the 
source does not easily permit separation of the locally 
generated field and the incident field for each cell. How- 
ever, such a separation is possible when is constructed 
by summation over a set of rays, which each start at a 
random position within the cell and point in a random 
direction. For ray i, call the incident radiation on the cell 
/o,i and the distance to the boundary of the cell dsi. The 
current level populations translate this distance into an 
opacity dr^, and give the source function Sui- The aver- 
age radiation field from N rays then follows from Eqs. ^ 
and (l7|), 



J. 



N 

^external 



J. 



local 



(14) 
(15) 



Here, Si, and dri contain both line and continuum terms, 
and /o,i includes the CMB. The radiation field is now 
the sum of the external (J^^t^nai) ^^^^ internal (4°"''') 
terms. Since the external term is evaluated using popula- 
tions from the previous Monte Carlo iteration (through 
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and Sui), this scheme is akin to accelerated A- iteration. 
Within Monte Carlo iterations, sub-iterations are used to 
find the locally self-consistent solution of Sui and for 

given jexternal^ 

The main computational cost of this strategy lies in fol- 
lowing a sufficient number of rays out of each cell through 
the source. Iteration on Eq. (|4|) is relatively cheap and 
speeds up convergence considerably in the presence of 
opaque cells. Fig. || illustrates this, by showing the evo- 
lution of the fractional error of the solution of the simple 
problem posed in §|| for optically thick HCO"*" and thin 
H^'^CO"'" excitation (for a fixed set of directions - see be- 
low). 

Population inversions require careful treatment in ra- 
diative transfer codes, since the associated opacity is neg- 
ative and the intensity grows exponentially. In general, an 
equilibrium will be reached where the increased radiation 
field quenches the maser. Since iterative schemes solve the 
radiative transfer before deriving a new solution for the ex- 
citation, the radiation field can grow too fast if population 
inversions are present. Our code handles negative opaci- 
ties by limiting the intensity to a fixed maximum which 
is much larger than any realistic field. Proper treatment 
requires that the grid is well chosen, so that masing re- 
gions are subdivided into small cells where the radiation 
field remains finite. Our code can deal with the small pop- 
ulation inversions that occur in many problems including 
the model presented in However, to m odel astrophysi- 
cal ma sers, specia lized codes are required ( ^paans and van 
Lange^ elde, 1992 , e.g.). 



3. 5. The role of variance in Monte Carlo calculations 

Because the Monte Carlo method estimates from a ran- 
domly chosen set of directions, the result has a variance, 
(T, which depends on the number N of included directions 
as CT cx 1/\/N. As explained above (§ |3.3| ), this variance 
is a strength rather than a weakness of the Monte Carlo 
method. Since it is not a priori known how many directions 
are required for a fiducial estimate of J^, this method auto- 
matically arrives at an appropriate sampling by increasing 
A'^ until the variance drops below a predefined value. 

The variance of a solution is usually estimated from 
the largest relative difference between subsequent itera- 
tions. In our implementation (see appendix), the number 
A^ of rays making up in a particular cell is doubled 
each time the variance in that cell exceeds a certain value; 
the variance is evaluated using the largest relative differ- 
ence between three subsequent solutions with the same A^. 
This cell-specific scheme ensures that the radiation field is 
sufficiently sampled everywhere, and at the same time pre- 
vents oversampling of cells which are close to LTE and/or 
weakly coupled to other regions. 

The variance as estimated from the difference between 
subsequent solutions only reflects the noise if random fluc- 
tuations dominate the difference. There will be systematic 



differences between subsequent solutions if these are still 
far from convergence. Therefore, many Monte Carlo meth- 
ods consist of two stages. In the first stage, a fixed number 
of photons will yield a roughly converged solution; in the 
second stage, the number of photons is increased until the 
noise becomes sufficiently small. 

In our implementation, this first stage consists of iter- 
ations with a fixed number of directions making up in 
each cell, Nq, which depends on the model. The directions 
are randomly distributed, but in each iteration, the same 
set of random directions is used by resetting the random 
number generator each iteration. Without random fluctu- 
ations in Jy, the difference between subsequent solutions 
purely reflects the evolution toward convergence. The flrst 
stage is considered converged when this 'noise' is a factor 
of ten smaller than the user-specifled level. 

For a sufficiently large Aq (typically a few hundred), 
the excitation in each cell now is close to the final solu- 
tion, except for imperfections in the sampling of Ji,. In the 
second stage, each iteration uses a different set of random 
directions to estimate Jj^: the random number generator is 
no longer reset. Based on the resulting variance, the num- 
ber of rays in each cell is doubled each iteration, until the 
noise on the level populations in each cell is below a given 
value. If A^o was initially insufficient, the variance will con- 
tain a significant contribution from systematic differences 
between iterations. Even though this will slow down the 
code by artificially increasing the number of rays in these 
cells as the code over-compensates the variance, ultimately 
the code will still converge to the correct solution. 

The separation between local and incident radiation 
fields in our method (§p.4[) keeps the system responsive to 
changes even in the presence of significant optical depth. 
This accelerates the convergence, but also increases the 
noise level. The literature mentions several methods to re- 
duce the noise of Monte Carlo methods, e.g., with a refer- 



ence field (Bernes, 1979; Ghoi ct al., 199£ ) or quasi- random 
instead of pseudo-random numbers ( luvela, 1997 ). These 
schemes are useful when assumptions about the solution 
are available, but may slow down convergence if the initial 
guess is far off. Since the 'first stage' described above and 
the presence of noise prevents Monte Carlo methods from 
'false convergence', we have not included any noise reduc- 
tion techniques in our code, to keep it as widely applicable 
as possible. 

3. 6. Implementation and performance characteristics 

Appendix ^ describes the structure of the program in de- 
tail, and provides a reference to a web site where we have 
made the source code of its spherically symmetric version 
publicly available. To test its performance. Appendix ^ 
compares results obtained with our code to those of other 
codes. 

The main part of the program deals with calculating 
the excitation through the source model. In a separate 
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Fi g. 4. Left: Densfty distribution (top) and temperature distribution (bottom) of the collapse model including rotation 
of §4.1. Middle: Resulting excitation temperature of HCO+ 1-0 and 3-2. The adopted grid is visible in these panels, 
with small cells at the center were the density is large, and larger cells in the outer regions of the object. Right: Images 
of integrated intensity of HCO+ 1-0 and 3-2 for an edge-on source orientation, after sampling on spatial frequencies 
corresponding to interferometric baselines between 15 and 300 m and subsequent image reconstruction. This results 
in synthesized beam sizes of 3" and 1" for the 1-0 and 3-2 lines, respectively, as indicated in the lower left corner of 
the panels. 



part, comparison to observations proceeds by integrating 
Eq. (|2|) on a grid of lines of sight for a source at a specified 
inclination angle and distance. The resulting maps of the 
sky brightness may be convolved with a beam pattern for 
comparison with single-dish data, or Fourier transformed 
for comparison with interferometer data. 

Appendix ^ describes tests of the validity of the re- 
sults of the program. We have also tested up to what 
optical depth the program can be used, and found that 
this depends on source model. These tests were done on a 
Sun Ultrasparc 2 computer with 256 Mb internal mem- 
ory and a 296 MHz processor. For a simple, homoge- 
neous HCO+ model with n = 10'^ cm~^ and T = 30 K, 
the code produces accurate results within an hour for 
values of 7V(HC0"'") up to ^ 10^'' cm~^, corresponding 
to T ~ 20,000 in the lowest four rotational transitions. 
Higher— J lines are less optically thick under these phys- 
ical conditions. For such opaque models, 'local' approxi- 
mations fail badly, because the excitation drops sharply 
at the edge of the model (Bernes 1979; Appendix]^. 

For a power-law, Shu-type model, performance is 
somewhat slower. The dense and warm region fills only 



a small volume, while its radiation has a significant in- 
fluence on the excitation further out, and modeling this 
effect requires a large number of rays. We have used the 
specific model from the Leiden workshop on radiative 
transfer (Appendix |^) for various values of the HCO+ 
abundance. Within a few hours, accurate results are ob- 
tained for values of HCO+/H2 up to 10^^, corresponding 
to r = 600 — 2000 in the lowest four lines. 



These test cases should bracket the range of one- 
dimensional models of interest. For two-dimensional mod- 
els, the limitations of present-day hardware are much more 
prohibitive. As a test case, we have used the flattened in- 
falling envelope model from I 
dances. Within 24 hours, 



4.1 for various HCO+ abun- 



the above machine provides a 
converged solution for HCO+/H2 up to 10~*, correspond- 
ing to a maximum optical depth of ~ 100. Realistic models 
often have higher opacities, and call for the use of paral- 
lel computers. However, as faster computers are rapidly 
becoming available, we expect that these limitations will 
become less relevant in the near future. For both one- and 
two-dimensional models, the second, ray-tracing part of 
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the code to create maps of the sky brightness takes only dered magnetic fields and/or rotation (Hogerheijde et al 



a fract on of the computer resources of the first part. 



f998, e.g.). These mechanisms probably influence the ac- 



3. 7. Alternative accelerators 

Another method to tackle slow convergence in the pres- 



ence of large opacities is core saturation (Rybicki, 1972; 
Doty and Neufeld, 1997), where photons in the optically 
thick line center are replaced by the source function and 
no longer followed, while photons in the still optically thin 
line wings which carry most information are more accu- 
rately followed. This scheme has been imple mented in a 
Monte Carlo program by Hartstein & Liseau ( 1998 ), but 
involves a choice where to separate the line core from the 
line wings. Since the effectiveness of the method depend 
on this choice, we have not included core saturation in our 
program. 

A completely different approach to accelerating con- 
vergence is to extrapolate the behaviour of the populations 
from the last f ew iterations. This so-called Ng accelera- 
tion ( Ng, 1974 ) is not implemented in our code, because 
extrapolating from an inherently noisy estimate may be 
dangerous. 

4. Astrophysically relevant examples 

A first example of the appl icat ion of our code has al- 
ready been given in §§|| and |3.4| . This model is a spheri- 
cally symmetric (one-dimensional) model; many authors 
have already illustrated the capability of Monte Carlo 
and other methods in solving one-dimensional problems 
( iBernes, 197"9| ; |Zhou, 1995| ; |Choi et al., 1995i e.g.). This 
section presents a number of astrophysically relevant ex- 
amples, again drawn from star formation studies, to illus- 
trate the distinguishing properties of our code - in ad- 
dition to accelerated convergence: the capability to cal- 
culate axisymmetric source models with a wide range of 
spatial scales, and the effects of dust continuum on radia- 
tive transfer and excitation. 

The models presented in the following sections all in- 
clude continuum radiation fields from dust. For these star- 
forming regions, we have chosen the model of Ossenkopf 
& Henning (1994) for the dust emissivity, which includes 
grain growth for a period of 10^ yr at an ambient density 
of 10^ cm~^ and thin ice mantles. Except for the calcu- 



lations in §4.3 where we specifically examine the effect of 



dust on the excitation including infrared transitions, only 
(sub) millimetre transitions were included in the excita- 
tion calculations which are not significantly influenced by 
the relatively weak continuum field. 

A young stellar object with rotation 

Observations of nearby young stellar objects often show 
flattened structures rath er tha n the spherical symmetry of 
models like that of Shu (1977), presumably caused by or- 



cretion behaviour, and may give rise to a rotationally sup- 
ported circumstellar disk. To test these ideas against ob- 
servations, cylindrically symmetric source models need to 
be considered. This section examines a model of protostel- 
lar collapse that includes flattening due to rotation follow- 



ing the description of Terebey, Shu, & Cassen (1984), and 
its appearance in aperture synthesis maps. 



The model of Terebey et al. ( 1984 ) treats rotat ion as 
a small perturbation to the solution of Shu (1977) for a 
collapsing envelope, joined smoothly to the description of 
a circumstellar disk by Cassen & Moosman (1981). In ad- 
dition to the sound speed and age, which we take identical 
to the values of §|| of a = 0.24 km s"^ and t = 1 x 10^ 
yr, this model is parameterized by a rotation rate Q. This 
gives rise to a centrifugal radius Rc = amQt^Q,^ /16, within 
which the infalling material accretes onto a thin disk. 
Here, toq = 0.975 is a numerical constant. We choose 
n = 5.9x 10-^3 s-\ so that Rc = 800 AU. We assume 
that inside Rc the material accretes onto a thin disk, and 
that most molecules rapidly freeze out onto dust grains 
(cf. §4.2). Effectively, we assume the region within Rc to 



be empty for this calculation. Fig. ^ (top left) shows the 
adopted density structure. All other parameters are simi- 
lar to the model of §||. 

To follow the power-law behaviour of the density in 
the model, a total of 15 x 15 cells are spaced exponen- 
tially in the R and z directions. To reach a final signal- 
to-noise ratio of 10, with 300 rays initially making up the 
radiation field in each of the cells, the Monte Carlo code 
requires 5 hours CPU time on a UltraSparc 10 to con- 
verge on the HCO"*" solution. For comparison, the opti- 
cally thin and more readily thermalized ^^CO problem 
takes only 10 minutes. The resulting excitation tempera- 
ture distribution (Fig. ^ middle panels) does not deviate 
much from that of the spherically symmetric model of §^ 
apart from the flattened distribution of the material at the 
center: rotation is only a small perturbation on the over- 
all source structure. As a result, the appearance is mostly 
unaffected in single-dish observations which do not resolve 
scales comparable to Rc where flattening is important. 

Millimetre interferometers, on the other hand, can re- 
solve these scales at the distances of the nearest star- 
forming regions, ~ 140 pc. Fig. ^ (right panels) shows 
the integrated emission in HCO+ J=l-0 and 3-2 after 
sampling at the same spatial scales as real interferometer 
observations and subsequent image reconstruction. Delays 
between 100-1000 ns were used, corresponding to angu- 
lar scales of 2"-40" for the 1-0 line and 0'.'6-13" for 3-2. 
Hence, emission on scales ^ 6000 AU (^ 2000 AU at 3-2) 
is filtered out. This results in a reconstructed ('cleaned') 
image that is dominated by the central, flattened regions 
of the envelope when the object is seen edge-on. Aperture 
synthesis observations of embedded protostars in Taurus 
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Fig. 5. The density, temperature, and molecular abundance distribution of the circumstellar disk model described 
in §4.2. Density and abundances are plotted on logarithmic scales; the temperature is plotted on a linear scale. The 



'superheated' layer of the disk model of Chiang & Goldreich is not included in our model because only a very small 
amount of dust and gas is present in this layer. Its 'backwarming' effect on the disk interior is included, however. 
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Fig. 6. The distribution of excitation temperatures in K for the J=l-0 and selected submillimetre lines of CO, 
HCO+, HCN, and CN. In all panels the grey scale ranges between K and 30 K on a linear scale. The adopted grid, 
exponentially spaced in radius and linearly in height, is reflected in the excitation distribution. 



show s imila r structures ( Ohashi et al., 1997 ; Hogerheijde 
et al, [99|, e.g.). 



4-2. A circumstellar disk 

Planetary systems form from the disks that surround 
many young stars dBeckwith, 1999| ; [Mannings etal.,2000| ) . 
Observational characterization of these disks is of prime 
importance to increase our understanding of the processes 
that shape planetary systems. Here, we present simu- 
lations of molecular line observations of a circumstellar 
disk around a T Tauri star as obtained with current and 
planned millimetre-interferometric facilities. 



The physical structure of the model disk is that of a 
passive accretion disk in vertical hydrostatic equilibrium 
as described by Chiang & Goldreich ( 1997 ). This descrip- 
tion includes the effect of 'backwarming' of the disk by 
thermal radiation of a thin, flared surface layer that in- 
tercepts the stellar light. The total amount of material in 
the superheated surface layer is too small to be detectable, 
but the overall effect of increased mid-plane temperature 
is significant. The effective temperature of the central star 
is 4000 K and its luminosity is 1.5 Lq. The outer radius 
of the disk is 700 AU, with a total mass of 0.04 Mq. 

The chemical structure of the disk follows Aikawa & 



Herbst (1999), who describe the radial and vertical com- 
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Fig. 7. Appearance of integrated intensity of selected lines in the circumstellar disk model at resolutions attainable 
with current and planned interferometric facilities. The first three columns from the left show CO 1-0 at 2" resoltuion, 
CO 1-0 at O'.'S resolution, and CO 6-5 at 0'.'2 resolution. A distance of 140 pc is adopted for the source, and inclination 
increases from 0° to 90° from top to bottom. The last three columns from the left show the emission in HCO+ 4-3, 
HCN 4-3, and CN 83/2 — > 2^/2 a-t 0'.'2 resolution. All panels are shown with a linear grey scale ranging from to the 
image maximum in K km s^^. Contours are drawn at 1%, 5%, 10%, 15%, 20%, 30%, 40%, 50%, 70%, and 90% of 
maximum. 



position of a flared disk. Freezing out of molecules onto 
dust grains is one of the dominant processes influencing 
the gas-phase composition in disks, and strongly depends 
on temperature and density. In the dense and cold mid- 
plane, many molecules will be depleted onto grains. How- 
ever, close to the star where temperatures are raised, and 
away from the midplane where densities are lower and de- 
pletion time scales longer, gas-phase abundances will be 
significant. In addition, ultraviolet radiation and X-rays 
may penetrate the upper layers of the disk, photodissoci- 
ating molecules and increasing the abundance of dissoci- 
ation products like CN. Fig. ^ shows the distribution of 
the density, temperature, and abundances of CO, IICO+, 
HCN and CN in the adopted model. We have used the 



results presented in Aikawa & Herbst (1999, their Figs. 6 



and 7; high ionization case) , and parameterized the abun- 
dances as function of scale height. 

For the Monte Carlo calculations, a gridding is adopted 
that follows the radial power-law density profile in 14 ex- 
ponentially distributed cells and the vertical Gaussian pro- 
file in 13 cells linearly distributed over 3 scale heights. 
Convergence to a signal-to-noise ratio of 10 requires ap- 
proximately 6 hours CPU time per model on an Ultra- 
Sparc 10 workstation, starting with 100 rays per cell and 
limiting the spatial dynamic range to 36 (i.e., the small- 



est cell measures 10 AU on the side). The resulting exci- 
tation and emission depends on the competing effects of 
increased abundance and decreased density with distance 
from the midplane. Fig. |6| shows the excitation tempera- 
ture of selected transitions and molecules. Fig. |^ shows a 
number of representative simulated observations, at res- 
olutions of 2", 0'.'5, and 0'.'2 attainable with current and 
planned (sub) millimetre interferometers. Van Zadelhoff et 
al. (in prep.) present a simpler analysis of similar models. 

4-. 3. A high mass young stellar object 

The above examples referred to the formation of stars 
with masses of up to few M0 and luminosities ^ 100 
L0. Stars of higher mass spend their first ~ lO'^ yr in 
envelopes of ~ 100 Mq. With their luminosities of 10"*- 
10^ L0, these stars heat significant parts of their en- 
velopes to several hundred K, shifting the peak of the 
Planck function to the wavelengths of the vibrational tran- 
sitions of many molecules. Stars of lower mass and lumi- 
nosity only heat small regions to a few hundred K, and 
the impact on the excitation is correspondingly smaller. 
As an example, Figure || shows two models of the HCN 
submillimctrc line emission, with and without pumping 
through the bending ('V2") mode at 14.02 /im. For com- 
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Fig. 8. Submillimetre emission lines in a 1" beam of HCN 
in the vibrational ground state (top panels) and first ex- 
cited state (bottom panels). Thin line: only coUisional ex- 
citation; thick line: model with pumping through the V2 
band at f 4 jim. 



putational convenience, only energy levels up to J—IQ 
within the first vibrationally excited and ground states 
have been included, including 1-type doubling in the ex- 
cited state. This doubling occurs due to the two possi- 
ble orientations of the rotational and vibrational motions 
with respect to each other. CoUisional rate coefficients be- 
tween rotational levels are from Green (f 994, priv. comm., 
see http://www.giss.nasa.gov/data/mcrateE); within vi- 
brational levels, they were set to 10"^'^ cm'^s"^ . The source 
model is that o f the y oung high-mass star GL 2136 by van 
der Tak et al. ( 2000| ). Based on its luminosity of 7 x 10"' 
L0 and dust mass of ~ 100 M©, the star has heated a 
region of radius ~ 3000 AU to ^ 100 K, making pumping 
through the 14 /zm HCN bending mode important. We 
have assumed that Tgas = Tdust, as is true for high density 



regions. The dust emissivity, from Ossenkopf & Henning 
( |1994| ), is the same as in the previous sections. As seen 
in Fig. H, the effect of dust is especially strong for the 
rotational lines within the V2=l band, which occur at fre- 
quencies slightly offset from the ground state transitions. 
These lines have indeed been detected towards similar ob- 
jects ( [Ziurys and Turnerl986 , e.g.). 

The shells around evolved stars is another area where 
inclusion of infrared pumping by dust is essential to 
understand 



the rotational line emission ( Ryde et al 
1999|7^.g.). Many molecules that are commonly observed 
through rotational lines at millimetre wavelengths have ro- 
vibrational bands in the mid-infrared, and can be pumped 
by warm dust. In a few cases, pumping through rotational 
lines at far-in frared wavelengths is im portant as well, for 
example CS ( Hauschildt et al., 1993|) and all hydrides. 



5. Conclusion 

This paper describes a computer code to solve the coupled 
problem of radiative transfer and molecular excitation for 
spherically and cylindrically symmetric source geometries. 
It is based on the Monte Carlo method, but incorporates 
elements from accelerated lambda iteration which greatly 
improve convergence in the presence of significant optical 
depth. In particular, the code separates excitation due to 
the local environment from the response to the radiation 
field after propagation through the source. This approach 
combines the flexibility of a Monte Carlo method with 
the reliability of accelerated lambda iteration. We expect 
our code to be a valuable tool in the interpretation of the 
wealth of data that current and future instruments op- 
erating from the millimetre to the infrared will yield. A 
number of examples (§^) already illustrates the applica- 
tions to problems in star formation studies. 
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Appendix A: The code 

Our code is applicable to a wide range of astrophysi- 
cal problems involving (sub) millimetre and infrared (^ 
10 /xm) spectral-line and continuum observations. The 
one-dimensional (spherically symmetrical) version of our 
code is publicly available for all interested researchers via 



http://astro.berkeley.edu/'^michiel. The two-dimensional 



(cylindrically symmetric) code is available on a collabora- 
tive basis through the authors (see the same website for 
contact information). This appendix gives a concise de- 
scription of the implementation of the accelerated Monte 
Carlo method. 



Fig. A.l gives an overview of the structure of our code. 



most notably water ( Hartstein and Liseau, 1998 ) 



which consists of two parts. The first runs the Monte Carlo 
simulation solving the radiative transfer and molecular 
excitation. The second part uses this solution to calcu- 
late the emission that would be observed from this source 
above the atmosphere and with perfect spatial and veloc- 
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Fig. A.l. Schematic outline of steps involved in our Monte 
Carlo calculations. 



ity resolution, given a source distance and, for cylindri- 
cally symmetric models, inclination. This latter part can 
also be used to calculate the continuum radiation emit- 
ted by the source. Its output format is that used by the 
MIRIAD package [Multichannel Image Reconstruction, 
Image Analysis, and Display; Sault et al. ( |1995| )]. This 
package, designed to analyse interferometer interferomet- 
ric spectral line data, includes many processing options 
such as convolution with a single-dish beam and mod- 
eling of aperture-synthesis visibilities, as well as a wide 
variety of imaging capabilities. MIRIAD also allows easy 
conversion to the ubiquitous FITS format and portability 
to other software packages. 

Both parts of the code are controlled by UNIX C-shcU 
scripts that extract information from the provided input 
and compile an executable code. In this way, the size of 
several arrays containing the source model, the coUisional 
rate coefficients, etc., can be adjusted to the required size, 
minimizing memory requirements. The source code is writ- 
ten in FORTRAN-77. 



Following the flow chart of Fig. A.l, the following steps 
describe the Monte Carlo part of the code in more detail. 



The code starts by reading a list of keywords, detailing 
the required signal-to-noise ratio on the level popula- 
tions, the initial number of photons in each cell (Nq), 
and pointers to the source model, the systematic veloc- 
ity field (if any), the description of the dust emissivity, 
and the molecular energy levels and collisional rate 
coefficients. The velocity field can be defined simply 
through the source model with each grid cell moving 
at a constant speed, or it can be a constantly varying 
function over each cell. The source model can be a se- 
ries of concentric shells covering a region from the ori- 
gin to a maximum radius, or a series of stacked cylin- 
ders fully covering a region out to a maximum radius 
and height. 

Collisional rate coefficients are available for many as- 
trophysically interesting species and common collision 
partners such as II2 in the J — and in the J — 1 
levels, e~, and He. Our code currently allows for two 
simultaneous collision partners, e.g., II2 and e~, each 
with its own density and temperature. For molecular 
ions such as 1100+, excitation due to collisions with 
electrons can be significant compared to collisions with 
H2 at fractional ionization 10~^). Often, listed rate 
coefficients are equivalent rates per H2 molecule in- 
cluding contributions from He at cosmic abundance. 
The results of our code, and any non-LTE calculation, 
sensitively depend on the quality of the rate coeffi- 



cients. Recently, Black (2000) discussed the need for 
good rate coefficients and the effects of other implicit 
assumptions of radiative transfer codes. 
In the first stage of the calculation, the radiation field 
is based on iVo rays per cell, each starting at a ran- 
dom position equally distributed over the cell volume, 
pointing in a random direction, and at a random fre- 
quency within 4.3 times the local line width around 
the local systematic velocity vector. The value of 4.3 
corresponds to the width where the line profile has 
dropped to less than 1% of its peak. In this stage, in 
each iteration the same series of random numbers is 
used, so that there are no random fiuctuations in the 
coverage of the radiation field. 

For each ray, the distance ds from the ray's origin to 
the nearest boundary of the cell along its randomly 
chosen direction is calculated. The incident radiation 
Iq along the ray then follows from integrating Eq. ^ 
in a stepwise manner from cell edge to cell edge, at- 
tenuating the contribution from each cell by all inter- 
vening cells, with the cosmic microwave background as 
a boundary condition. The only quantity that changes 
when stepping through a cell is the direction, and pos- 
sibly the magnitude, of the systematic velocity vector, 
which enters Eq. ^ through the line profile function 
0(1^). Changes of 0(i^) within cells are tracked by sub- 
dividing the integration into small steps as needed. 
Armed with the set of {ds , , Iq) for each ray, the 
radiation field in the cell follows from Eq. (|4[). A 
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consistent solution of this equation and the level pop- 
ulations (Eq. |l^) quickly follows from iteration to a 
relative accuracy of 10~® in the populations. Limita- 



tions on masering are discussed in Sect. 3.4. 
The first stage of the code repeats items 2-4 until the 
largest relative fractional difference between the pop- 
ulations in all cells of three subsequent solutions is 
ten times better than ultimately required. Since the 
angular sampling is the same in each iteration, these 
differences are free of random noise but might not ad- 
equately sample all directions and frequencies. 
The second stage of the code proceeds along similar 
lines as the first stage, but with a different set of ran- 
dom numbers in each iteration. The only other differ- 
ence is, that each time the maximum fractional error 
in the populations in a cell exceeds the requested ac- 
curacy, the number of rays in that cell is doubled. This 
stage lasts until all cells comply with the required ac- 
curacy, after which the solutions are written out to a 
file. 





CO 2-1 






□ Bernes (1979) 
This work 





Radius (m) 




Velocity (km s ^) 



Fig. B.l. Excitation temperature of CO J = 1 ^ and 
2 — > 1 as a function of radius, and integrated line profiles. 
Open symbols are results by Bernes (1979), the sofid fines 
are our calculations. 



The second part of the program calculates the emis- 
sion distribution on the sky for a given source distance 
and inclination by simple ray tracing. The output from 
the Monte Carlo code forms the input for this ray-tracing 
code. Since it uses much of the same code as the Monte 
Carlo part, geometry and radiative transfer being the 
same, it is not further discussed here. 

Appendix B: Comparison with other codes 

This section describes two cases to test our code 
against well-documented calculations with Monte Carlo 
codes from the literature. For further tests, we re- 
fer the reader to the web-page collecting a num- 
ber of standard test cases, which has resulted from 
the 1999 workshop on Radiative Transfer in Molecu- 
lar Lines at the Lorentz Center of Leiden University 
(http://www.strw.lcidcnuniv.nl/~radtran6). 



B.l. Bernes' CO cloud 

In his seminal paper on Monte Carlo methods for radiative 
transfer and molecular excitation, Bernes ( 1979| ) presents 
a constant-density, constant-temperature, optically thick 
cloud model. The density of the cloud, = 2000 cm~^, 
is below the critical density of the CO transitions, and the 
excitation is dominated by radiative trapping. The excita- 
tion temperatures of the CO transitions drop off rapidly 
in the outer regions of the cloud. This necessitates fine 



sampling of these regions. Figure B.l shows that our code 
reproduces the original results within the accuracy of our 
and Bernes' calculations. This simple model forms a criti- 
cal test for the code's ability to correctly handle excitation 
by radiative trapping. The total run time for the model 
was approximately 5 minutes on a UltraSparc 10 worksta- 



— Our code. ID 

Our code, 2D 

□ Choi et al. (1995) 




-2.5 -2 -1.5 -1 



-2.5 -2 -1.5 -1 
log(radius/pc) 



-2.5 -2 -1.5 -1 



Fig. B.2. Excitation temperature of selected CS and 
H2CO lines versus radius. Open symbols are results by 
Choi et al. (1995), the thin solid lines are the results from 
our spherically symmetric code. The thick solid lines show 
the results from our cylindrically symmetric code. 



tion, using the same coUisional rate coefficients as Bernes 
(Table 0). 

B.2. Model for B 335 by Choi et al. 

Another critical element of any radiative transfer code is 
its ability to correctly deal with systematic velocity fields. 
The inside-out collapse model as outlined in §^ is well 
suited for such a test, because of its wide range in ve- 
locities from zero to many times the turbulent line width 
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combined with significant optical depth. As a test case, 
we calculate the populations and the emergent spectrum 
of several CS and H2CO lines, following the model for 
the young stellar object B 335 of Choi et al. (1995). This 
model is similar to that of § |^, with a = 0.23 km s^^ and 
t = 1.3 X 10*5 yr. The turbulent line width is 0.12 km s^^. 
Only the temperature structure is different from § ^ Choi 
et al. use continuum observations to constrain the temper- 
ature distribution, which we follow as closely as possible 
from their Fig. 3. 



Fig. B.2 compares the resulting excitation temperature 



distribution with the results of Choi et al. The agreement 
is very good for CS, where we used the same collisional 
rate coefficients (Table |l|). For H2CO the agreement is 
less favourable, but we were unable to use the exact same 
rate coefficients. Simple tests show that the deviation is 
comparable to what can be expected from the difference in 
the molecular data. This variation corresponds to a 10% 
difference in the emergent line profiles. 

Fig. B.2 also plots the excitation temperatures ob- 



tained for the same model but using the cylindrically sym- 
metric code. Both codes clearly give consistent answers; 
the small 'wiggles' in the excitation temperatures as func- 
tion of radius in the output of the cylindrically symmetric 
calculation can be attributed to geometrical defects when 
trying to fit a sphere in a series of stacked cylinders. 
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